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ABSTRACT 

The  rationale  and  numerical  technique  of  embedding  an  oceanic 
bulk  mixed  layer  model  with  a  multi-level  primitive  equation  model 
is  presented.  In  addition  to  the  usual  prognostic  variables  that 
exist  in  a  multi-level  primitive  equation  model,  the  embedded 
model  predicts  the  depth  of  the  well  mixed  layer  as  well  as  the 
jumps  in  temperature  and  velocity  that  occur  at  the  base  of  that 
layer.  The  depth  of  the  mixed  layer  need  not  coincide  with  any  of 
the  fixed  model  levels  used  in  the  primitive  equations  calcula- 
tions. 

In  addition  to  advective  changes,  the  mixed  layer  can  deepen 
by  entrainment  and  it  can  reform  at  a  shallower  depth  in  the  ab- 
sence of  entrainment.  When  the  mixed  layer  reforms  at  a  shallower 
depth,  the  vertical  profile  of  temperature  below  the  new,  shal- 
lower mixed  layer  is  adjusted  to  fit  the  fixed- level  structure 
used  in  the  primitive  equations  calculations  using  a  method  which 
conserves  heat,  momentum  and  potential  energy.  Finally,  a  dynamic 
stability  condition,  which  includes  a  consideration  of  both  the 
vertical  current  shear  and  the  vertical  temperature  gradient,  is 
introduced  in  place  of  the  traditional  "convective  adjustment". 

A  two-dimensional  version  of  the  model  is  used  to  test  the  em- 
bedded model  formulations  and  to  study  the  response  of  the  ocean 
to  a  stationary  axisymmetric  hurricane.  The  model  results  indi- 
cate a  strong  interdependence  between  vertical  turbulent  mixing 
and  advection  of  heat. 


1.   INTRODUCTION 

The  value  of  numerical  models  for  the  study  of  ocean  dynamics 
has  become  widely  recognized  in  recent  years  (National  Academy  of 
Sciences,  1975;  Kraus,  1977;  WMO/ICSU,  1977).  However,  no  single 
model  has  been  developed  that  can  resolve  the  full  range  of  space 
and  time  scales  commonly  observed  in  the  ocean.  Thus,  the  charac- 
teristics of  each  numerical  model  differ  considerably  and  they  are 
usually  tailored  to  the  space  and  time  scales  of  the  oceanic  phe- 
nomena being  studied.  The  time  scales  to  which  various  models  have 
been  applied  range  from  decades  for  the  deep  ocean  circulation  to 
seasonal,  synoptic  and  diurnal  for  near-surface  phenomena. 

Numerical  models  of  the  ocean  can  be  separated  into  two  broad 
categories.  These  are  the  basin-wide  or  global  ocean  models  which 
attempt  to  simulate  the  ocean  circulation  on  seasonal  or  longer 
time  scales,  and  the  one-dimensional  models  which  attempt  to  simu- 
late the  locally  forced  response  at  seasonal  or  shorter  time 
scales.  The  former  type  of  model  is  referred  to  as  an  ocean  gen- 
eral circulation  model  (OGCM)  and  the  latter  type,  which  usually 
focuses  on  vertical  mixing  processes,  is  referred  to  as  a  mixed 
layer  model  (MLM) .  The  most  recent  comprehensive  reviews  of  OGCM's 
and  MLM's  are  those  by  Holland  (1979),  Robinson,  Harrison  and 
Haidvogel  (1979),  Bryan  (1979),  Haney  (1979)  and  Garwood  (1979). 

The  present  study,  which  addresses  the  problem  of  embedding  a 
one-dimensional  mixed  layer  model  into  a  numerical  model  of  the 
general  circulation,  was  motivated  by  several  factors.  First,  the 
solution  of  this  problem  was  recently  set  as  an  important  goal  by 
an  international  group  of  scientists  studying  ocean  models  and 


their  relation  to  climate  (WMO/ICSU,  1977).  In  addition,  there  are 
a  wide  variety  of  important  ocean  problems  involving  intermediate 
and  mesoscale  motions  which  could  be  addressed  with  a  realistical- 
ly embedded  mixed-layer  ocean  circulation  model.  These  problems 
include  studies  of  open  ocean  fronts,  coastal  and  equatorial  up- 
welling  and  the  ocean's  response  to  intense  atmospheric  storms.  It 
is  also  anticipated  that  an  embedded  mixed  layer  ocean-circulation 
model  will  eventually  have  a  number  of  important  applications  in 
the  problem  of  operational  ocean  prediction. 

In  this  paper,  the  rationale  and  numerical  technique  of  embed- 
ding a  model  of  the  mixed  layer  into  an  ocean  circulation  model  is 
presented.  The  resulting  embedded  model  has  been  tested  in  a 
variety  of  ocean  studies.  The  simulation  of  the  ocean's  response 
to  hurricane  forcing  is  presented  here  because  the  large  horizon- 
tal variations  in  advection  and  mixing  in  such  a  case  (Madala  and 
Piacsek,  1977)  provide  a  severe  test  of  the  embedded  model.  Model 
simulations  of  the  oceanic  response  to  less  extreme  forcing,  which 
illustrate  more  subtle  aspects  of  the  embedding,  will  be  published 
later.  Our  purpose  is  not  to  intercompare  possible  mixing  para- 
meterizations,  but  rather  to  describe  the  embedding  of  a  particu- 
lar bulk  mixed  layer  model  within  a  multi-level  ocean  circulation 
model . 
2.   MODEL  EQUATIONS  AND  BOUNDARY  CONDITIONS 

In  this  section,  the  equations  and  boundary  conditions  for  the 
embedded  mixed-layer  ocean  circulation  model,  as  applied  to  the 
hurricane  forcing  problem,  are  presented.  Details  of  the  numerical 
formulation,  however,  are  given  in  the  Appendix. 


The  ocean  is  taken  to  be  hydrostatic  and  incompressible  and 
the  density  is  a  linear  function  of  temperature  alone.  For  compu- 
tational economy  and  because  hurricane  forcing  is  nearly  axisymme- 
tric,  the  ocean's  response  is  assumed  to  be  axisymmetric.  In  addi- 
tion, the  Coriolis  parameter  is  assumed  to  be  constant.  The 
governing  equations  are  written  in  radial  coordinates  as 
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where  u,   v  and  w  are  the  radial,    tangential  and  vertical  compo- 
nents of  velocity,    respectively;  p  is  pressure;  T  is  temperature; 
p   is  density;   A^.  and  A„  are  coefficients  of  horizontal  eddy  vis- 
cosity and  conductivity,    respectively;    a  is  the  coefficient  of 


thermal  expansion;  and  p  is  the  density  of  water  at  the  reference 
temperature,  T  .  The  values  of  constants  used  in  the  hurricane 


forcing  experiments  are  given  in  Table  1.  The  terms  (w'  u' ) , 
(w' v' )  and  (w'T')  represent  ensemble  average  turbulent  fluxes  and 
are  described  below.  In  (2.1)- (2.3),  the  advection  terms  have 
been  written  in  flux  form  using  the  continuity  equation  (2.4). 

The  above  model  equations  are  applied  to  an  oceanic  region  ex- 
tending from  the  surface  to  a  depth  of  H  =  964  m  and  from  the 
storm  center  (r  =  0)  to  a  radial  distance  of  R  =  450  km.  There 
are  no  fluxes  of  mass,  momentum  or  heat  across  the  boundaries  ex- 
cept at  the  surface  where  a  steady  wind  stress  and  upward  heat 
flux  are  prescribed.  The  specific  conditions  are 
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(w'T')  =  Q(r)/pQC 
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The  forcing  functions  Q,  T  and  t  are  shown  in  Fig.  1.  The 

r      9 

stress  is  calculated  from  a  prescribed  wind,  the  magnitude  of 
which  increases  linearly  with  r  from  the  eye  wall  of  the  storm 

(r  =  4.5  km)  to  the  radius  of  maximum  winds  (r  =  45  km);  de- 

-1/2 
creases  like  r     from  there  to  360  km;  and  then  decreases 

linearly  to  zero  at  the  boundary  (450  km) .  This  wind  profile 

results  in  the  wind  stress  components  shown  in  Fig.  1  and  a 

curl,  3(t  r)/rar,  which  is  zero  inside  the  eye  of  the  storm; 
9 

positive  and  increasing  linearly  from  r  =  4.5  km  to  45  km;  zero 
from  45  km  to  360  km  and  negative  from  r  =  360  km  to  the  bound- 
ary. The  maximum  tangential  stress  of  nearly  36  dPa  corresponds  to 
a  maximum  surface  wind  speed  of  50  m  s  .  The  small  inward  radial 
stress  component,  t  <  0,  results  from  a  cross  isobar ic  surface 
inflow  angle  of  20  degrees.  The  effect  of  other  choices  of  wind 
stress  profiles  will  be  reported  elsewhere. 

Putting  w  =  0  at  the  top  of  the  ocean,  as  in  (2.8),  is  known 
as  the  rigid-lid  approximation  and  it  requires  the  divergence  of 
the  vertically  averaged  motion  to  be  zero.  In  the  present  model, 
we  simplify  the  dynamics  even  further  by  assuming  that  the  ver- 
tically averaged  motion  itself  is  zero.  This  simplification  is 
justified  in  the  present  case  because  the  shear  currents  spun  up 
by  the  surface  stress  are  much  greater  than  the  dynamically  in- 
duced barotropic  currents.  Subtracting  the  vertical  average  of 
(2.1)  and  (2.2)  from  (2.1)  and  (2.2)  respectively  gives  the  pre- 
diction equations  for  the  vertical  shear  currents 
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where   (     )    denotes  a  vertical  average  over  the  total  ocean  depth 
H,    p  is  p-p,   and    (u,v)    are  now  the  shear  currents  having   zero  ver- 
tical average.   By  integrating    (2.5)    and  using    (2.6),    p  is  obtained 
as  a  function  of  T  alone,   thus 


Po(l^t(T-To))gdA 


I/./.. 

-H  z 


(l-a(T-T  ))gdX]dz    .      (2.13) 


The  turbulence  closure  method  described  below  requires  a  prog- 
nostic equation  for  the  depth  of  the  well-mixed  layer,  denoted  by 
h.  This  equation  is  derived  by  integrating  the  continuity  equation 
(2.4)  over  the  mixed  layer,  and  applying  the  boundary  condition 
w  =  0  at  z  =  0.  Thus, 
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Denoting  an  average  over  h  by  <(  )>  =  tr  /  (  )  dz  and  noting  that 


dh   8h  ,     _8h 
dt  ""  dt   U-h  8r  ' 


(2.14)  reduces  to 
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The  right  hand  side  of  (2.15),  denoted  w  ,  is  an  entrainment  ve- 
locity which  is  related  to  the  net  acquisition  of  mass  in  the 
mixed  layer.  Its  parameterization  is  described  below.  The  six 
equations  (2.11)f  (2.12),  (2.3),  (2.13),  (2.4)  and  (2.15),  along 
with  the  boundary  conditions  (2. 7) -(2. 10)  represent  a  closed  sys- 
tem of  equations  for  the  six  dependent  variables  u,v,T,p,w  and  h 


provided  the  turbulence  quantities,  (w' u'),   w  etc.,  are  para- 
meterized. The  values  of  the  turbulent  fluxes  of  heat  and  momentum 


at  the  top  and  bottom  of  the  model  ocean  are  given  by  the  boundary 
conditions  (2.9)  and  (2.10).  The  parameterization  of  the  turbulent 
fluxes  within  the  water  column  is  given  in  the  following  sections. 
3.   VERTICAL  DIFFUSION 

The  classical  solution  to  the  problem  of  parameterizing  the 
vertical  turbulent  fluxes  of  heat  and  momentum  is  the  posing  of 
eddy  viscosity  and  eddy  conductivity  formulations  such  that 


9u 
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w'T'  -  -  K^ 


This  is  reasonable  if  gradient  diffusion  is  physically  plausible, 
and  if  the  values  of  K,  and  K_  are  known  to  be  approximately  con- 
stant. This  approximation  has  been  shown  to  produce  realistic 
vertical  profiles  of  temperature  in  the  upper  ocean  below  the 
depth  of  immediate  influence  of  the  stirring  action  of  the  atmo- 
spheric forcing.  Therefore,  this  parameterization  is  employed 
at  all  depths  below  the  ocean  surface  mixed  layer,  with  KM  and  PL, 
both  equal  to  0.5  cm  s  .  It  is  possible  that  the  dynamic  insta- 
bility parameterization  described  later  could  alone  provide  the 
desired  mixing  below  the  turbulent  boundary  layer. 

The  above  parameterization  of  the  vertical  fluxes  is  not  sa- 
tisfactory, however,  between  the  ocean  surface  and  the  top  of  the 
stable  thermocline.  In  this  region  of  intense  turbulent  mixing, 
the  appropriate  values  of  PL.  and  K„  would  be  two  to  three  orders 
of  magnitude  larger  and  time  dependent.   Even  if  realistic  values 


10 


(3.1) 


of  K^  and  K  can  be  computed  from  scaling  arguments  or  from  more 
sophisticated  turbulence  closure  assumptions,  the  most  important 
problem  is  to  determine  the  depth  below  which  the  intensity  of  the 
turbulence  becomes  insignificant.  The  position  of  this  disconti- 
nuity is  at  z  =  -h,  the  depth  of  the  bottom  of  the  ocean  surface 
turbulent  boundary  layer  or  mixed  layer.  Thus  an  entrainment 
model  is  needed  to  compute  h(r,t)  . 
4.   ENTRAINMENT  AND  MIXED  LAYER  MODEL  FORMULATION 

In  the  Garwood  (1977)  model,  the  buoyancy  flux  at  the  base  of 
the  turbulent  boundary  layer  is  modeled  according  to 

ag  GPY7)   .  =  -  <w'w'>1/2  <  E  >/h  .  (AtL) 

-n 

This  expression  is  derived  from  the  local  (z  =  -h)  turbulent 
kinetic  energy  budget.   In  the  entrainment  zone,  the  buoyant  damp- 


ing attributable  to  (w'T1)  ,  is  assumed  to  be  balanced  by  the  con- 
vergence of  turbulent  kinetic  energy  flux.  This  convergence  term 
is  modeled  as  a  function  of  the  turbulent  kinetic  energy  compo- 
nents and  the  vertical  distance  over  which  this  energy  must  be 


transported.  Here  <w*w' >,  and  <E>  are  the  bulk  (mixed  layer  aver- 
age) values  of  the  vertical  component  and  total  turbulent  kinetic 
energies  respectively.  The  prognostic  equations  for  these  mixed 
layer  turbulence  variables  are  derived  using  the  bulk  second  order 
closure  methods  of  Garwood  (1977)  and  are 

3   „  .=.  ,ns  3 


(h<E>/2)  =  mu/  -  ctgh  (w'T')_h/(2Ri  ) 


-  (<E>1/2  +  fh)  <E>  (4.2) 
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|^  (h<v'w'>/2)    =  agh    ((w'T')_h  -    (w'T')Q)/2 

+  (<E>  -3<wTwr>)    <E>1/2  -   (<E>1/2  +  fh)    <E>/3    ,        (4.3) 

1/2 
In   (4.2)    and    (4.3),    u*  =    (t/p  )  with  x  the  magnitude  of  the 

* 

surface  stress  and  p  the  (constant)  air  density,  and  Ri  =  aghAT/ 

a 

2     2 
(au  +  av  )  is  the  bulk  Richardson  number.  The  quantities  AT,  au 

and  av  are  the  temperature  and  velocity  jumps,  or  differences, 

between  the  values  in  the  mixed  layer  (assumed  to  be  vertically 

uniform)  and  the  underlying  stable  water.  These  bulk  turbulent 

kinetic  energy  equations  are  solved  algebraically  by  assuming  a 

quasi-steady  state  for  the  turbulent  kinetic  energy  budgets: 


Ij-  (h  <E>)  -  0  ;  -^  (h  <w^7">)  -  0 


This  assumption  causes  the  computed  values  of  <w'w' >  and  <E>  to  be 
in  phase  with  the  surface  boundary  conditions — the  wind  stress, 


2 


p  u*  ,  and  the  effective  surface  buoyancy  flux,  ag(w'T')  .  In  a 
study  of  nonstationarity  including  a  diurnal  cycle  (Garwood  and 
Yun,  1979),  this  approximation  was  found  to  be  sufficiently  accu- 
rate for  time  scales  greater  than  a  few  hours.  The  motivation  for 
this  quasi-steady  state  assumption  is  a  substantial  gain  in  numer- 
ical efficiency.  The  time  step  is  limited  more  by  the  requirement 
of  resolving  the  diurnal  forcing  than  for  accuracy  in  solving  the 
turbulent  kinetic  energy  equation.  This  assumption  also  makes  the 
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boundary  layer  model  more  compatible  with  the  numerical  require- 
ments of  the  ocean  circulation  model. 

For  some  applications,  the  bulk  Richardson  number  term  in 

* 
(4.2)  may  also  be  neglected  if  Ri  »  1.  However,  this  term 

should  be  retained  in  situations  in  which  strong  shears  might  be 

expected  across  the  entrainment  zone.  Niiler  (1975)  and  DeSzoeke 

and  Rhines  (1976)  demonstrate  that  this  term  may  be  significant  at 

the  outset  of  a  deepening  event,  provided  the  initial  mixed  layer 

depth  is  small. 

Once  the  entrainment  buoyancy  flux  is  known  from  (4.1),  the 

downward  fluxes  of  momentum  associated  with  entrainment  at  the 


base  of  the  mixed  layer,  (w'u1),  and  (w'v')_h/  are  given  by 

-  (^T)_h  =  weAu  =  -  (wTTT).h  |r  (4.1a) 

"  Cw^T)   =  wAv  =  -  GTT)^  §  ,  (4.1b) 


where  the  entrainment  velocity, 


-h 
e      Ai 


(4.1c) 


The  above  mixed  layer  formulation  was  selected  for  coupling 
with  an  ocean  circulation  model  because  of  its  general  applica- 
bility and  because  it  has  been  extensively  tested  and  calibrated 
at  several  weather  ship  stations  (Garwood  and  Adamec,  1980) .  The 
general  applicability  of  the  above  model  is  attributed  to  two 
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important  features:  (1)  An  annual-period  cyclical  steady  state  is 
possible  because  of  a  dissipative  limitation  to  deepening  that  is 
dependent  upon  planetary  rotation;  (2)  Entrainment  flux  is  a  func- 
tion of  implicitly  computed  values  of  both  horizontal  and  vertical 
components  of  turbulent  kinetic  energy.  As  a  result,  the  fraction 
of  turbulent  kinetic  energy  that  is  available  for  mixing  at  the 
base  of  the  mixed  layer  is  not  a  constant  but  is  dependent  upon 
both  the  surface  buoyancy  flux  and  the  surface  friction  velocity. 
These  desirable  properties  are  obtained  using  second  order  closure 
techniques  without  sacrificing  the  bulk  model's  numerical  effi- 
ciency which  is  necessary  for  practical  application  of  the  final 
coupled  system  with  large  two  and  three  dimensional  arrays. 
Boundary  layer  deepening  by  entrainment 

The  entrainment  heat  flux,  (w'T')  ,  ,  is  determined  algebra- 
ically from  (4.1),  and  the  steady  state  forms  of  (4.2)  and  (4.3) 


using  given  values  of  u^,  (w'T')  ,  h  and  Ri  .  This  heat  flux  is 
then  imposed  upon  the  given  temperature  profile.  The  new  h  and 
the  new  momentum  and  temperature  profiles  are  unique  solutions 
provided:  (i)  the  mixed  layer  (0  >  z  >  -h)  is  homogeneous;  (ii) 
the  value  of  h  is  made  just  large  enough  to  prevent  unstable  den- 
sity profiles  (9T/az  >  0);  and  (iii)  heat  is  conserved.  Note  that 
the  depth  of  the  bottom  of  the  mixed  layer  is  not  forced  to  coin- 
cide with  any  of  the  prescribed  model  levels  or  interfaces  between 
levels.  As  a  result,  there  is  added  vertical  resolution  in  the 
vicinity  of  the  mixed  layer  interface.  This  extra  resolution  is 
very  desirable  because  the  vertical  structure  of  the  thermocline 
determines  the  potential  energy  of  the  upper  ocean  which  is 
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important  for  both  the  dynamics  of  the  boundary  layer  and  the 
larger  scale  baroclinic  circulation. 
Boundary  layer  shallowing 

Equally  important  in  a  mixed  layer  model  of  general  applica- 
bility is  the  prediction  of  boundary  layer  shallowing.  This  phe- 
nomenon occurs  when  there  is  inadequate  vertical  turbulent  energy 
to  transport  heat  down  to  the  base  of  the  existing  mixed  layer. 


The  (new)  depth  at  which  the  vertical  flux  (w'T')  vanishes  es- 
tablishes the  bottom  of  the  turbulent  boundary  layer  and  what  will 
be  the  base  of  the  new  mixed  layer.  In  the  present  model,  this 


occurs  when  <w' w'>  approaches  zero.  Then  (4.1)  no  longer  applies 
and  the  steady  state  forms  of  (4.2)  and  (4.3)  reduce  to  the 
following  diagnostic  equations  for  the  new  mixed  layer  depth 

0  =  mu^3  -  (<E>1/2  +  fh)  <E>  (A. 2a) 

0  =  -agh(wTTr)o/2  +  <E>3/2  -  (<S>1/2  +  fh)  <S>/3  .  (4.3a) 

Without  the  term  involving  the  planetary  rotation,  or  if  the  down- 
ward surface  buoyancy  flux  is  dominant,  h  is  proportional  to  the 


Obukhov  length  scale,  L  =  u*  /(ag(w'T')  ).  On  the  other  hand,  if 


(w'T1)  =  0,  h  is  proportional  to  u^/f.  In  general,  however,  the 
depth  of  the  shallowing  mixed  layer  is  a  function  of  the  two  non- 
dimensional  parameters,  h/L  and  hf/u*. 

In  a  numerical  model,  the  formation  of  a  new,  shallower  mixed 
layer  is  potentially  a  problem  because  the  interfacial  structure 
of  previously  created  mixed  layers  can  not  be  exactly  preserved 
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with  a  relatively  coarse  vertical  grid.  In  the  present  model, 
this  potential  problem  is  avoided  using  a  numerical  procedure  (see 
Appendix)  that  conserves  heat,  momentum  and  potential  energy, 
while  minimizing  changes  in  the  vertical  density  structure. 
Although  potential  energy  is  conserved  during  the  shallowing  pro- 
cess, no  attempt  is  made  to  exactly  conserve  mean  kinetic  energy. 

5.  DYNAMIC  STABILITY  CONDITION 

In  the  derivation  of  the  above  bulk  mixed  layer  model,  it  is 
assumed  that  the  turbulent  boundary  layer  is  dynamically  unstable 
and  that  the  underlying  water  column  is  dynamically  stable  as 
measured  by  the  gradient  Richardson  number,  Ri.  In  the  entrap- 
ment zone  (-h  <_  z  _<  -h-d) ,  dynamic  stability  is  assumed  to  be  near 
neutral;  that  is  Ri  =  Ri   =  1/4.  This  condition  is  insufficient 
in  itself  to  determine  the  entrainment  fluxes  or  the  rate  of  mixed 

layer  deepening;  however,  it  does  lead  to  a  measure  of  the  thick- 

2     2 
ness  of  the  zone:  d  =  Ri~D  Uu  +  av  )agAT.  Deepening  of  the  mixed 

CR 

layer  depends  upon  the  added  intermittent  vertical  flux  of 
turbulent  kinetic  energy  from  above. 

As  mentioned  earlier,  dynamic  instability  can  sometimes  occur 
below  the  mixed  layer,  and  this  will  enhance  the  vertical  trans- 
port of  both  heat  and  momentum.  Therefore,  in  the  present  model, 
dynamic  stability  is  enforced  throughout  the  water  column  by  im- 
posing those  vertical  fluxes  of  heat  and  momentum  between  model 
levels  which  are  necessary  to  prevent  the  gradient  Richardson  num- 
ber from  going  below  the  critical  value,  i.e.  Ri(z)  >  Ri,—,.  This 

—   CR 

"dynamical  adjustment"  is  a  simple  generalization  of  the  "convec- 
tive  adjustment"  parameterization  which  is  commonly  used  in  most 
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ocean  general  circulation  models.  Details  are  given  in  the 
Append  ix . 

It  should  be  pointed  out  here  that  the  above  dynamical  sta- 
bility condition  is  different  from  that  of  Pollard,  et  al.  (1973) 

and  Thompson  (1976) .  In  that  model,  deepening  is  dependent  upon 

* 

stability  as  determined  by  the  bulk  Richardson  number,  Ri  .  This 

is  equivalent  to  the  present  requirement  that  there  be  a  constant 
gradient  Richardson  number  in  the  entrainment  zone  only  if  the 
thickness  of  the  entrainment  zone  divided  by  the  mixed  layer 
depth,  i.e.  d/h,  is  a  constant.  There  is  no  physical  justifica- 
tion for  such  an  assumption. 
6.   METHOD  OF  EMBEDDING 

The  system  of  equations  for  the  embedded  model  includes  a 
prognostic  determination  of  the  mixed  layer  depth.  The  mixed 
layer  is  defined  here  as  being  the  upper  layer  of  the  ocean  where 
turbulent  mixing  is  dominant.  The  temperature  and  horizontal 
velocities  are  required  to  be  independent  of  depth  in  this  region, 

One  of  the  more  important  parameters  in  predicting  the  evolu- 
tion of  the  mixed  layer  is  the  dynamic  stability  at  its  base. 
Since  density  is  assumed  to  be  a  function  of  temperature  alone, 
the  temperature  jump  at  the  base  of  the  mixed  layer  is  critical. 
The  dynamic  portion  of  the  model  is  a  level  model  that  predicts 
the  average  of  a  quantity  between  two  specified  depths.  The  base 
of  the  mixed  layer,  however,  will  rarely  lie  exactly  on  a  model 
level.  Requiring  a  sufficiently  fine  vertical  resolution  to  pre- 
vent errors  in  the  thermal  and  potential  energy  balances  due  to 
the  truncation  of  the  mixed  layer  depth  to  a  model  level  is 
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unreasonable  both  in  terms  of  computer  storage  and  computational 
time.  The  problem  addressed  here  is  to  provide  for  a  meaningful 
communication  between  the  fixed-level ,  dynamical  part  of  the 
model,  and  the  mixed  layer  part  of  the  model. 

The  coupling  of  the  dynamical  and  mixing  processes  in  the 
model  may  be  considered  in  two  stages.  First,  the  changes  in  the 
upper  ocean  due  to  advective  and  diffusive  processes  are  calcu- 
lated in  the  dynamical  part  of  the  model  and  put  into  appropriate 
form  for  the  mixed  layer  model.  Second,  the  changes  due  to  sur- 
face fluxes  and  entrainment  mixing  are  calculated  in  the  mixing 
part  of  the  model  and  transmitted  to  the  dynamical  part  of  the 
model.  In  both  of  these  stages,  a  special  treatment  is  required 
for  the  model  level  which  contains  the  base  of  the  mixed  layer. 

At  the  start  of  a  given  timestep,  changes  due  to  advective  and 
diffusive  processes  (including  the  dynamic  adjustment  described  in 
Section  5)  are  calculated  for  all  the  GCM  layers  in  the  dynamical 
part  of  the  model  and  the  prognostic  variables  (u,v,T  and  h)  are 
stepped  forward  in  time  (see  the  Appendix  for  details  of  the  time 
differencing  scheme) .  The  special  treatment  for  the  layer  con- 
taining the  new  mixed  layer  depth  is  now  described  with  the  aid  of 
Fig.  2  in  which  the  dashed  line  is  the  vertical  structure  used  by 
the  dynamical  part  of  the  model  and  the  solid  line  is  the  vertical 
structure  used  by  the  mixing  part  of  the  model.  Both  the  dynami- 
cal and  the  mixing  part  of  the  model  use  the  same  structure  at  all 
levels  which  lie  entirely  within  the  mixed  layer.  The  variable  e, 
represents  u,  v,  or  T.   If  the  base  of  the  mixed  layer  lies  within 
the  k  th  model  level,  the  value  of  5  in  the  mixed  layer  portion  is 
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denoted  by  5  ;  the  value  just  below  the  mixed  layer  is  denoted  by 
5  •  the  integrated  average  over  the  model  level  is  denoted  by  \; 
while  the  mixed  layer  depth  is  denoted  by  h.  The  change  55  due 
to  advective  and  diffusive  processes  is  assigned  the  value  calcu- 
lated for  the  GCM  layer  immediately  above  (i.e.  the  layer  bounded 
by  Z,  _  .     and  Z,       .     in  Fig.  2).  The  base  of  the  mixed  layer  is 
constrained  to  always  lie  at  or  below  the  base  of  the  first  model 
level.  The  associated  change  6£  is  then  calculated  from  the  re- 
quirement that  the  weighted  average  of  6?1  and  6E,     be  equal  to  the 
change  6£  predicted  for  the  layer  by  the  dynamical  part  of  the 
model .  Thus , 


AZk  65  =  21   6C1  +  Z2  6?2  ,  (5.1) 

where  Z,  is  the  portion  of  the  layer  above  the  new  mixed  layer  and 
Z.  +  Z9  =  AZ,  .  In  the  circulation  model,  the  advective  and  diffu- 
sive changes  are  based  on  the  average  properties  in  each  layer. 
However,  there  can  be  a  large  velocity  and  temperature  jump  across 
the  interface  h.  The  use  of  a  layer  mean  advecting  velocity  and 
layer  mean  temperature  in  the  dynamical  part  of  the  model  neces- 
sarily introduces  an  error  in  the  calculation  of  the  change  &l, 
which,  by  our  formulation,  is  introduced  into  the  change  assigned 
to  6£„.  The  error  in  6£  will  be  smallest  when  the  depth  h  is 
close  to  one  of  the  interfaces  between  GCM  layers,  and  it  will  be 
largest  when  h  is  near  the  center  of  a  GCM  layer.  An  analogous 
but  slightly  more  complex  formulation  is  used  whenever  the  dynam- 
ically predicted  sh  causes  the  base  of  the  mixed  layer  to  move 
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into  an  adjacent  model  layer.  Before  proceeding  to  the  mixing 
portion  of  the  timestep,  the  values  of  e;  throughout  the  mixed 
layer  (including  5,)  are  replaced  by  the  vertical  average  of  5 
over  the  new  mixed  layer  depth.  In  the  particular  hurricane  simu- 
lation shown  below,  the  small  amount  of  turbulent  kinetic  energy 
that  would  be  expended  (or  released)  in  performing  this  vertical 
mixing  has  been  neglected  in  the  turbulent  kinetic  energy  balance. 
This  is  clearly  justified  in  the  present  case  since  the  amount  of 
turbulent  kinetic  energy  generated  by  shear  stress  production  is 
so  large. 

The  second  stage  of  the  coupling  is  the  mixing  portion  of  the 
timestep.  During  this  stage,  the  case  with  entrainment  mixing  oc- 
curring in  association  with  mixed  layer  deepening  must  be  distin- 
guished from  the  case  in  which  the  layer  reforms  at  a  shallower 
depth.  Consider  first  the  entrainment  mode  in  which  the  mixed  layer 
model  (Section  4)     predicts  an  entrainment  heat  flux,  (w'T1 )_, 
<  0.  The  mixed  layer  deepening,  6h  >  0,  during  the  timestep  At  is 
calculated  from  (4.1c) .  In  this  case,  a  volume  of  water  (per  unit 
area)  equal  to  <5h  and  having  property  5_  is  entrained  into  the 
mixed  layer  having  volume  (per  unit  area)  equal  to  h  and  property 
£, .  This  results  in  a  new  mixed  layer  of  depth  h  +  6h  and  property 
E,     +  6£,,  where 


fir   .  -1_^.     .  (P. 2) 


Entrainment  mixing  produces  no  change  in  £  below  the  new  mixed 
layer  depth.   If  5,  >  £„,  which  is  always  true  for  the  temperature 
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in  the  model,  then  6S.  <^  0.  If  the  entrainment  deepening  does  not 
cause  the  base  of  the  mixed  layer  to  move  into  an  adjacent  model 
layer,  then  the  magnitude  of  the  property  jump  across  the  mixed 
layer,  | £,  -  59l,  will  be  decreased.  In  the  numerical  model,  a 
separate  array  of  K?   values  must  be  stored  to  permit  the  jump 
values,  5,  -  5_,  to  be  always  specified  for  the  mixing  process. 
This  value,  and  the  mixed  layer  depth  h  are  the  only  additional 
variables  which  must  be  retained  in  coupling  the  mixing  model  and 
the  dynamical  model  (see  Appendix  for  details) . 

Consider  now  the  case  of  mixed  layer  shallowing.  In  this  case 
the  diagnostic  relations  (4.2a)  and  (4.3a)  are  solved  for  the  new 
mixed  layer  depth.  The  heat  and  momentum  fluxed  across  the  sur- 
face during  the  timestep  are  distributed  over  the  new,  shallower 
depth  to  produce  a  new  value  of  5, .  A  new  £,<-,   value  is  determined 
from  the  prior  5, .  The  next  advective  timestep  operates  with  the 
new,  shallower  mixed  layer.  Some  information  is  lost  regarding 
the  deeper  structure  because  provision  is  made  for  but  one  active 
mixed  layer.  To  preserve  some  of  the  information  regarding  the 
density  jump  at  the  base  of  a  prior  mixed  layer,  the  vertical 
density  profile  is  adjusted  to  conserve  both  heat  and  potential 
energy  whenever  the  mixed  layer  reforms  at  a  shallower  depth  (see 
Appendix  for  details) . 
7.   APPLICATION  WITH  STATIONARY  HURRICANE  FORCING 

The  purpose  of  this  section  is  to  demonstrate  the  capability 
of  the  embedded  mixed  layer-ocean  circulation  model  to  simulate 
realistic  upper  ocean  changes  in  response  to  hurricane  forcing. 
Such  strong  forcing  as  that  of  a  hurricane  is  used  because  a  large 
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horizontal  gradient  in  the  advective  and  mixing  response  is  expected 
and  this  will  provide  an  extreme  test  of  the  model.  No  attempt  is  made 
to  compare  the  results  of  the  present  embedded  model  to  the  results  of 
other  models  utilizing  different  mixing  parameteri zations .  Such  a  com- 
parison is  beyond  the  scope  of  the  present  paper. 

The  classical  model  simulation  of  the  ocean's  response  to  hurricane 
forcing  is  that  of  O'Brien  and  Reid  (1967)  using  an  axisymmetric,  one- 
layer  model.  The  extension  of  this  type  of  model  to  two  horizontal 
dimensions  is  given  by  Chang  and  Anthes  (1978) ,  while  the  response  of  a 
multi-level  circulation  model  is  described  by  Madala  and  piacsek 
(1977) .  As  a  result  of  these  and  other  studies,  certain  features  of 
the  dynamical  response  of  the  ocean  to  a  stationary  hurricane  are  well 
known.  The  currents  in  the  upper  layers  of  the  ocean  tend  to  circulate 
in  the  same  sense  as  the  hurricane  winds.  Although  the  hurricane  wind 
field  is  directed  not  only  cyclonically  about  the  storm  but  also 
radially  inward,  the  net  radial  movement  of  surface  layer  water  is  out- 
ward due  to  the  Coriolis  effect.  A  vertical- radial  circulation,  con- 
sisting of  intense  upwelling  within  the  radius  of  maximum  hurricane 
winds  and  weaker  downwelling  at  greater  radial  distance   is  produced 
by  mass  continuity.  As  a  result  of  this  circulation,  the  ocean  tem- 
peratures near  the  center  of  the  storm  decrease  throughout  the  stable 
part  of  the  water  column  that  lies  below  the  mixed  layer.  The  ocean 
surface  temperatures  are  also  affected  if  the  upwelling  water  reaches 
the  surface. 

There  is  also  an  important  mixing  effect  due  to  the  hurricane  force 
winds,  as  indicated  by  the  model  studies  of  Elsberry,  Fraim  and 
Trapnell  (1976)  and  Chang  and  Anthes  (1978)  .  Mechanical  and  convective 
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generation  of  turbulent  kinetic  energy  in  the  upper  ocean  tend  to 
achieve  maximum  values  at  the  radius  of  maximum  winds.  As  a  result, 
the  relatively  warm  surface  water  is  cooled  and  the  coler  water  imme- 
diately below  the  initial  mixed  layer  depth  is  warmed  as  it  is  en- 
trained into  and  is  mixed  with  the  surface  water.  The  combination  of 
the  advective  effects  and  mixing  effects  described  above  is  clearly 
quite  complex  and  requires  high  vertical  and  horizontal  resolution  to 
simulate. 

To  simplify  the  testing  of  the  present  model,  the  forcing  due  to  a 
stationary,  axisymmetric  hurricane  is  used.  A  4.5  km  horizontal  grid 
spacing  is  utilized  in  order  to  have  10  grid  points  within  the  radius 
of  maximum  winds  which  is  located  at  45  km  (Fig.  1)  .  The  ocean  is 
initially  at  rest  with  a  uniform  surface  temperature  of  30  C  and  a 
uniform  mixed  layer  depth  of  60  m.  To  reduce  the  initial  shock,  the 
hurricane  forcing  is  slowly  spun  up  to  full  intensity  over  the  first 
12  h. 

The  evolution  predicted  at  12  h  and  24  h  with  the  embedded  mixed 
layer-ocean  circulation  model  is  presented  in  Figs.  3a-d  and  Figs. 
4a-d,  respectively.  Because  of  the  explicit  mixed  layer,  the  surface 
inputs  of  heat  and  momentum  are  distributed  over  a  deep  layer.  Thus, 
the  maximum  v  component  at  12  h  (Fig.  3a)  is  rather  small  considering 
the  large  magnitude  of  the  surface  stress  near  the  center  (Fig.  1)  .  At 
12  h  there  is  a  deep  layer  of  cyclonically  rotating  water  in  the  mixed 
layer  and  a  very  small  and  weak  anticyclonic  circulation  below  the 
mixed  layer.  A  region  of  relatively  weak  outward  flow  is  predicted 
(Fig.  3b)  for  the  mixed  layer  in  the  inner  region.  The  compensating 
inward  flow  below  is  a  maximum  near  the  top  of  the  thermocline.  As  is 
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indicated  by  the  packing  of  the  w-isolines  at  the  center  of  the  hurri- 
cane in  Fig.  3c,  the  associated  vertical  motion  at  12  h  is  concentrated 
in  a  very  narrow  zone.  The  temperature  field  at  12  h  in  Fig.  3d  is 
beginning  to  show  the  effect  of  the  strong  upwelling  near  the  center. 
Toward  the  center,  the  isotherms  are  tilted  upward  at  progressively 
larger  angles.  The  29  C  isotherm  has  been  lifted  about  15  m  near  the 
storm  center. 

At  24  h  the  tangential  (v)  components  within  the  mixed  layer  near 
the  center  are  only  slightly  stronger  than  at  12  h.  A  maximum  value  of 
less  than  1  m  s   is  qualitatively  in  agreement  with  the  largest  cur- 
rents measured  as  hurricane  Eloise  passed  an  instrumented  oil  rig  off 
the  coast  of  Louisiana  (Forristall,  1980).  Although  the  strength  of 
the  hurricane  forcing  used  here  is  somewhat  larger  than  in  Eloise, 
there  is  a  compensating  effect  due  to  the  larger  layer  depth  in  the 
model  compared  to  the  preEloise  ocean  conditions  of  between  30  and  45 
m.  At  larger  r,  there  is  a  region  of  weak,  anticyclonic  flow  that  is 
produced  as  the  water  moves  outward  while  tending  to  conserve  angular 
momentum.  The  radial  extent  of  the  outward  flow  near  the  surface  in 
Fig.  4b  is  much  greater  than  at  12  h.  The  maximum  of  this  radial  com- 
ponent is  nearly  as  large  as  the  corresponding  tangential  component,  so 
that  the  current  direction  is  approximately  40°  relative  to  a  circle, 
or  60  to  the  right  of  the  surface  stress.  The  compensating  inflow  oc- 
curs over  a  deep  layer,  although  the  maximum  is  clearly  located  in  the 
upper  layers  of  the  thermocline  beneath  the  radius  with  maximum  stress. 
The  primary  upwelling  at  24  h  remains  at  the  center  of  the  hurricane 
and  extends  to  great  depths  (only  the  upper  layers  of  the  model  are 
shown  here)  .  A  much  weaker  descent  occurs  at  large  r  to  compensate  for 
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the  surface  flow,  which  has  reached  the  outer  boundary  by  24  h.  The 
primary  feature  in  the  temperature  field  at  24  h  (Fig.  4d)  is  the  ex- 
pansion of  the  zone  of  upwelled  water,  which  has  surfaced  near  the 
center  and  has  been  advected  outward. 

The  advective  effect  near  the  center  of  the  hurricane  is  so  domi- 
nant that  it  is  difficult  to  isolate  the  mixing  effect.  Deviations  of 
the  predicted  temperature  fields  at  12  h  and  24  h  from  the  initial 
values  are  shown  in  Figs.  5a  and  5b,  respectively.  After  12  h,  the 
maximum  temperature  decrease  occurs  near  the  storm  center  in  associa- 
tion with  the  strong  upwelling  there,  while  a  region  of  temperature 
increase  occurs  near  the  outer  boundary  due  to  downwelling.  There  is  a 
second  zone  of  increased  temperatures  at  the  base  of  the  mixed  layer, 
which  has  deepened  as  a  result  of  both  entrainment  and  vertical  advec- 
tion.  A  similarly  shaped  region  of  temperature  increase  is  found  some- 
what farther  from  the  center  of  the  24  h  field  shown  in  Fig.  5b.  Such 
a  region  of  temperature  decrease  above  the  initial  mixed  layer  base  and 
increase  below  it  is  an  indication  of  entrainment  mixing.  It  is  clear 
from  the  above  temperature  changes  that  in  the  present  case  of  station- 
ary forcing,  advective  effects  dominate  the  heat  budget.  If  the  hurri- 
cane was  translating,  upwelling  at  a  point  would  be  followed  by  down- 
welling  and  the  advective  effects  would  have  a  tendency  to  cancel.  In 
such  a  case,  the  mixing  effects  would  be  relatively  more  important  in 
producing  a  permanent  change  in  the  thermal  structure. 

The  explicit  treatment  of  the  dynamics  of  the  mixed  layer  in  the 
embedded  model  results  in  a  temperature  field  that  is  similar  to  obser- 
vations of  oceanic  response  to  hurricanes  (e.g.  Leipper,  1967).  The 
purpose  of  the  limited  test  shown  here  is  only  to  demonstrate  that  the 
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embedding  procedure  proposed  here  appears  to  be  capable  of  treating 
both  the  mixing  effects  and  the  advective  effects  in  a  numerically  ef- 
ficient and  physically  realistic  manner.  No  inference  should  be  made 
regarding  the  viability  of  different  mixing  parameter izat ions. 
8.   CONCLUDING  REMARKS 

This  paper  has  addressed  the  problem  of  developing  a  dynamical 
ocean  circulation  model  capable  of  treating  both  advective  and  mixing 
effects  in  response  to  realistic  atmospheric  forcing.  The  solution 
presented  here  is  the  embedding  of  a  bulk  mixed  layer  model  within  a 
multi-level  primitive  equation  model,  in  addition,  the  traditional 
"convective  adjustment"  often  used  in  such  primitive  equation  models 
has  been  replaced  by  a  more  general  dynamical  stability  condition.  The 
resulting  model  is  an  attempt  to  combine  the  advantages  of  layer  models 
and  level  models  into  a  single  model. 

Other  approaches  to  the  problem  are  clearly  possible.  One  alterna- 
tive approach  is  to  use  a  higher  order  diffusion  model  in  which  the 
vertical  eddy  viscosity  and  conductivity  are  expressed  as  functions  of 
the  model  variables.  This  approach  has  been  utilized  by  Madala  and 
piacsek  (1977)  using  the  popular  formulas  of  Munk  and  Anderson  (1948)  . 
With  such  an  approach,  high  vertical  resolution  is  necessary  in  order 
to  represent  the  sharp  transition  between  the  highly  turbulent  surface 
layer  and  the  weakly  turbulent,  highly  stratified  region  below. 
Another  alternative  approach  is  to  develop  a  multi-layer  model  which 
includes  turbulent  exchange  between  the  layers.  A  bulk  mixed  layer 
could  be  considered  as  the  uppermost  layer  in  such  a  model.  One  of  the 
results  of  the  present  study  suggests  that,  with  this  approach,  there 
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should  also  be  a  turbulent  transfer  of  momentum  between  the  layers,  an 
effect  that  has  not  been  generally  included  in  layer  models  before. 

The  example  of  the  oceanic  response  to  a  stationary,  axisymmetric 
hurricane  examined  here  is  an  extreme  test  of  the  present  embedded 
model  because  both  the  advective  and  the  mixing  effects  are  very 
strong.  Additional  studies  with  the  present  model  are  currently  in 
progress.  These  include  the  ocean's  response  to  a  variety  of  strong 
atmospheric  forcing  situations  and  the  response  of  an  upper  ocean 
density  front  to  local  atmospheric  forcing.  A  three  dimensional  ver- 
sion of  the  model  is  being  used  to  study  a  variety  of  oceanic  meso- 
scale  features  as  well  as  large-scale  interannual  variability  (anoma- 
lies) in  the  central  North  pacific  Ocean.  These  studies  suggest  that 
the  embedded  model  is  capable  of  simulating  the  ocean' s  response  to  a 
wide  range  of  oceanic  and  atmospheric  forcing  conditions.  Other  tests, 
including  comparisons  with  observational  data,  are  planned  to  further 
validate  the  present  model.  As  a  result,  it  now  appears  that  the  model 
can  serve  as  a  useful  tool  to  increase  our  understanding  of,  and 
eventually  our  ability  to  predict,  a  variety  of  oceanic  phenomena  that 
involve  both  advective  and  mixing  processes. 
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APPENDIX 

Space  finite  differencing 

The  numerical  scheme  makes  use  of  centered  space  differencing 
on  a  staggered  grid.  In  the  notation  which  follows,  the  letters  j 
and  k  are  integers.  The  vertical  index  is  k  and  the  variables  u, 
v,  T  and  p  are  defined  at  level  k,  k=l,  ...K,  while  the  vertical 
velocity  w  is  defined  at  the  K-l  interfaces  between  the  levels  and 
also  at  the  top  and  bottom  of  the  ocean  where  it  is  zero. 
Standard  difference  notation  is  utilized  in  which  for  example 

y*i+i/2,j =  «i+i,j  -  ii,j and  rqX)i+i/2,j =  1/2  (qi+i,j + 
qi,j'- 

Let  AZ.  ,  k=l,...K,  be  the  specified  thickness  of  the  k  th 
level  and 


Z  ,  =  -  Y     AZ   ;  k=l,  ...  K  (A3) 


be  the  value  of  Z  at  the  base  of  the  k  th  level  (Z,/?  =  0  and 
Z„  , /?  =  -H) .  Average  layer  thicknesses  are  defined  by 


AZ    -   CA.Z\,_,  ;  k-l,  ...  K-l  ,  (A2) 

and  average  layer  heights  are  defined  by 

Z,.  -    (//},..  ;    k-l,  ...  K  .  (A3) 
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The  finite-difference  analog  of  (2.13)  is 

Pk  =  Pk  "  P  .  (A4) 

where  the  wavy  overbar  denotes  the  finite-difference  vertical 
average 

K 


('">  4E  (  >kAzk  •  (A5) 


k-l 


and 

k-1 

Thus,  p,  is  the  thermally  induced  hydrostatic  pressure  field  which 
has  a  zero  vertical  average. 

Let  Ar  be  the  (constant)  horizontal  grid  size.  Then  the  veloc- 
ity components  u  and  v  are  defined  at  the  center  of  the  grid  in- 
tervals located  at  the  points 

r  =  \   Art-  (j-1)  LT    ;  j-1,  ...  J  .  (A7) 

The  variables  T,  p  and  w  are  defined  at  the  J-1  interfaces  and  at 

both  ends  of  the  domain  (r=0  and  R) .  These  points  are  located  at 

r,  -  0 

r,fi,  -  (rr),H,  :        j-1,  ...  J-I  \  •  (A3) 

r J*  "  R 
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In  addition,  the  space  differencing  scheme  makes  use  of  the  fol- 
lowing definitions  of  the  grid  volumes. 


ir^,k  =  2  rl  Ar  AZk 


7rj+^,k  "  Li+h  "^   ~k 


r_,.  Ar  AZ,.  ;  j-1,  ...  J 
j=l,  ...  J 


*J.k"  ^J.k  : 


(A9) 


is 


The  finite-difference  analog  of  the  continuity  equation  (2.4) 


6r(u'V;>k  +  6Z(w*>j-R5,k  =  °  ;   J=0'  ■"  J  ■ 


(A10) 


where 


i .  ,    =  u .  ,  r  .  AZ, 
J,k      j,k  j    k 


Vs.w*  =  Vf.w*  r**  Ar 


(All  ) 


The  finite  difference  analog  of  the  boundary  conditions  which  are 
needed  to  apply  (A10) ,  and  the  equations  to  follow,  on  the  bound- 
aries are  given  in  the  next  section. 

The  finite-difference  analog  of  the  thermal  energy  equation 
(2.3)  and  the  prognostic  equation  for  the  mixed  layer  depth  (2.15) 
are 
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3t    (Tj4^,k   irj+^,k) 


(A12) 


8t   (  V>  rj4-  Ar> 


-   ^(hu)*^,  +  (we)j+1.  t^h  Ar 


(A13) 


where 


(T4l:  ^^j.k^k 


Oto) 


-,z  * 


(Tr) 


Am  fl    (T) .   .    r .    KL  /Ar 


(A1A) 


V><     ±    ;,-i;     1,11 


(hu) 


(^H ..„iUb  Vi  Ar 


(Jir>     u^ 


The  finite-difference  analog  of  the  radial  component  of  the 
equation  of  motion  (2.11)  is 


7) 


-,-,:   (u.  ,  tt.  ,  ) 

j  »  K  J  >  '^ 


+  <f  +  VJ.l:/rj)vj.k,'j.fc*  '5r®j,krj  AVpc 
+  6r(ur)  ,,k  "  V^j.k  '  W»jik 


(A15) 
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where, 


/-r. 


AZ,    K 
k=l 

V.  .   11.  ,        (t  )  .  IT.  , 

p  H 


The  finite-difference  analog  of  the  tangential  component  of  the 
equation  of  motion  (2.12)  has  a  similar  form. 
Boundary  conditions 

The  lateral  boundary  conditions  will  only  be  given  for  r=0 
since  the  conditions  at  r=R  are  similar.  The  variables  and  T,p 
and  w  are  defined  at  grid  points  which  coincide  with  the  lateral 
boundaries  while  (u,v)  are  defined  at  points  located  one-half  grid 
distance  inside  the  boundary.  Thus,  the  variables  r,  ,»,  T,  ,~  .  , 
p...-  .  and  w  ,  ,   ,  are  located  at  the  origin  while  r. ,  u..  .  and 
v,  .  are  located  one-half  grid  distance  from  the  origin. 

The  zero  normal  flux  condition  at  r  =  0  is  maintained  by  im- 
posing the  following  antisymmetric  condition  on  the  normal  flux: 


u0,k  =  "  ul,k  '  (A17) 
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(A16) 


while  the  condition  of  no  lateral  eddy  flux  of  heat  and  momentum 
at  r  =  0,   which  is  used  in   (A12)    and    (A15) ,    is  given  by 


<Tr)o,k  "  -  <Tr>l,k 


^ 


(ur) 


i,k 


(vr) 


3,k 


(M8) 


Equations  (A17)  and  (A18)  are  the  finite-difference  analogs  of 
(2.7). 

The  finite-difference  analogs  of  the  boundary  conditions  on  w 
in  (2.8)  are 


<  ,,K^ 


0  . 


The  finite-difference  analog  of  the  boundary  conditions  on  the 
surface  fluxes  in  (2.9)  and  (2.10)  are 


(AW) 


(w'T')   .  ,  =  (()   ,./p  C)  r.,,,  Ar 


(w'vr).  ^   -  ((\0,/Pj  r  Ar 


(A20) 


and 


=  ( 


J  A  '"'2 


(A21) 


respectively,  where  the  forcing  functions  Q,  t  and  x  are  taken 
from  Fig.  1. 
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Vertical  diffusion 

In  this  section,  and  the  ones  to  follow,  the  numerical  methods 
used  to  calculate  the  turbulent  fluxes  in  (A12) ,  (A13)  and  (A15) 
are  given.  The  vertical  flux  of  heat  and  momentum  due  to  eddy 
diffusion  is  calculated  at  the  top  and  bottom  (interface)  of  the 
model  level  which  contains  the  mixed  layer  depth  and  at  the  top 
and  bottom  of  every  model  level  which  lies  entirely  below  the 
mixed  layer.  That  part  of  the  turbulent  fluxes  in  (A12)  and  (A15) 
which  is  due  to  vertical  diffusion  is  calculated  from 


(w"i")j^|k+li  -  -  Kj  y~)j+;„W,  rj+!j  Ar/AZ^ 


.'- 


rajM.  --W^.k^4*"^ 


with  a  similar  equation  for  (w*  v' )• 
Boundary  layer  entrainment 

This  section  concerns  the  numerical  treatment  of  entrainment 
by  the  embedded  mixed- layer  model.  In  a  time  step  At,  the  amount 
of  mixed  layer  deepening  and  the  associated  layer  temperature 
change  depend  upon  both  the  prior  thermal  structure  and  the  en- 
trainment heat  flux  (w'T')_h,  computed  by  solving  (4.1)  simul- 
taneously with  (4.2)  and  (4.3).  Following  (6.2),  heat  is  verti- 
cally redistributed,  and  the  changes  in  mixed  layer  temperature 
and  depth  are  found  as  illustrated  by  Fig.  Al  in  which 

i)    6T  =  (wT"")^  At/h 
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(A22) 


ii)  &h   is  determined  by  requiring 
area  B  =  area  A  =  h  6T  . 
Only  after  applying  this  entrainment  mixing  is  the  new  mixed  layer 
temperature  further  adjusted  to  account  for  the  net  surface  heat 


flux,  (w'T')   (although  (w*T')  helped  to  determine  the  value  of 
5T  above  through  its  appearance  in  the  turbulent  kinetic  energy 
equation,  (4.3)).  Once  the  new  mixed  layer  depth  is  determined, 
the  horizontal  components  of  mixed  layer  velocity  are  adjusted  so 
as  to  conserve  momentum  and  preserve  vertical  homogeneity  in  the 
mixed  layer. 
Boundary  layer  shallowing 

When  the  depth  of  the  turbulent  boundary  layer  decreases,  it 
leaves  behind  a  remnant  mixed  layer  structure.  In  a  numerical  si- 
mulation of  this  event,  unless  the  base  of  the  older  mixed  layer 
exactly  coincides  with  an  interface  between  two  model  levels,  some 
information  must  be  lost  because  only  one  mixed  layer  depth  is  re- 
tained. This  loss  of  information  is  a  potentially  serious  problem 
because  such  remnant  mixed  layer  structures  influence  the  poten- 
tial energy  in  the  column.  In  the  present  model,  when  the  mixed 
layer  depth  is  predicted  to  decrease,  from  h  to  h'  say,  the  tem- 
perature in  the  model  levels  which  lie  partly  or  entirely  within 
the  mixed  layer  are  adjusted  along  with  the  temperature  in  the 
model  level  immediately  below  the  mixed  layer  depth  (Fig.  A2) . 
This  adjustment  is  given  by 
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1  Vfc=l,k 


k+1 


Tk-,1  +  " 


h(Tk  "  Vi)[J\  - h  +  V 


k     n; 


Tk+i  + 


m 


k 


-  Max(h'f    D        ) 


(A23) 


(A24) 


(A25) 


v^ere  k  is  the  model  level  which  contained  the  original  mixed 
layer  depth  before  shallowing,   T.    and  h  are  the  model   temperatures 
and  mixed  layer  depth  before  shallowing,   while  T.  '    is  the  adjusted 
model   temperatures  after  shallowing.      It  can  be  shown  that  the 
above  adjustment  of  temperatures  conserves  both  heat  and  potential 
energy. 

The  procedure  is  illustrated  in  Fig.  A2.     A  slight  vertical 
redistribution  of  heat  enables  the  conservation  of  potential 
energy.     In  the  case  when  the  new  turbulent  boundary  layer,  h1 ,    is 
shallower  than 

k-1 

Vi  -  Z  Az*  > 

£=1 

the  new  mixed  layer  depth  will  not  be  apparent  in  the  temperature 
profile  until  the  turbulent  boundary  layer  is  sufficiently  warmed 
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by  a  net  flux  of  heat  and  radiation  at  the  surface.  In  the  illus- 
tration, area  C  +  area  D  =  area  E.  In  this  purposely  exaggerated 
example,  the  surface  temperature  is  seen  to  rise  a  slight  amount. 
However,  over  a  cycle  of  deepening  shallowing  and  deepening  again, 
there  will  be  no  net  change  in  surface  temperature. 
Dynamic  stability  condition 

This  stability  condition,  which  is  a  simple  generalization  of 
the  familiar  "convective  adjustment",  is  imposed  by  examining  the 
value  of  the  Richardson  number,  Ri,  defined  at  the  interface  be- 
tween model  levels  given  below.  If  the  inequality  in  (A26)  is  not 
satisfied,  the  temperature  and  the  velocity  components  in  the  two 
adjacent  levels  are  "dynamically  adjusted"  in  accordance  with  the 
following  four  assumptions: 

(i)    heat  is  conserved, 

(ii)   momentum  is  conserved, 

(iii)   dynamic  stability  is  imposed  between  the  two  levels  by 
requiring  that  the  Richardson  number 


ctg  Az  ,,  AT 

Ri -^ j  >    (Ri>   -  .25  ,  (A26) 

(Au)   +  (Av)   "" 

(iv)   the  mixing  ratios  for  heat  and  momentum  are  equal,  i.e. 
\ _k  ==  Xk-H     k-tl  ^   \ _k  = 

a  UH-1  "  "k+1  a  \1  VK   ,  ^J±T.5±l  ,  .  (A27) 

uk+l  - "       vk  •  ■ ;     vk.i  - "' 
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where  the  prime  denotes  the  new,   adjusted  value.     In  these  expressions, 
AT,   Au  and  Av  are  the  vertical  differences  of  heat  and  momentum, 
T.    -  T.        etc.,   while  T,  u  and  v  represent  the  weighted  vertical 
average  of  the  two  levels,   0.5(TkAZk  +  \+ih\+,) /h\+1/2  etc*     The 
above  stability  condition  is  calculated  at  the  velocity  points  in  the 
model,   and  therefore,  because  of  the  staggered  grid,   the  temperatures 
appearing   in   (A26)    and    (A27)    are  two  grid-point  averages.     Conse- 
quently,  a  weighted  average  of  the  temperature  adjustment  determined 
from  the  above  conditions   is  used  to  compute  the  adjustments  that  are 
applied  at  the  temperature  grid  points.     With  some  straightforward 
algebraic  manipulations   it  can  be  shown  that  the  unique  solution  to 
conditions   (i)    through   (iv)    is 


T'        =  T  +     ybT 
k 


K+i " f "  *7AT 


u,'       =  u  +     yAu 


lk+l 

k 


u  -  CyAu 


v  +     YAv 


(A28) 


where 


V,'  ■-■    V    -    T/y-Av 

k+J 


Ri 


Y   p 


"(Ri).CR     ^,) 


and 
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If  the  mixing  ratios  in  (A27)  are  denoted  by  m  ,  it  can  easily  be 
shown  that 

m  =  1  ;  Ri  <  0  .  (A29) 


Complete  mixing  between  two  levels,  as  defined  by  m=l,  occurs 
only  if  the  value  of  the  Richardson  number  before  mixing  is  less 
than  or  equal  to  zero. 
Time  integration 

Time  integration  is  performed  using  the  leapfrog  scheme  with 
an  Euler-backward  scheme  introduced  at  the  start  and  at  every  11 
timesteps  thereafter  to  remove  any  solution  separation  in  time. 
The  forward  differencing  scheme,  however,  is  used  for  the  eddy 
diffusion  and  turbulent  mixing  terms.  As  indicated  in  Section  6, 
the  calculation  of  the  advection  and  diffusion  terms  is  split  from 
the  calculation  of  the  turbulent  mixing  terms. 

Suppose  the  dependent  variables  u,v,T  and  h  are  known  at  the 
two  consecutive  time  levels  n-1  and  n.  When  using  the  leapfrog 

scheme,  the  first  step  is  to  calculate  the  pressure  p  and  the  ver- 

* 
tical  velocity  w  using  (A4)  and  (A10)  and  the  values  of  T,  u  and 

v  at  time  level  n.  Next,  the  dynamical  terms  on  the  right  hand 
side  of  (A12) ,  (A13)  and  (A15)  are  calculated,  also  using  the 
variables  at  time  level  n.  These  dynamical  terms  include  the  hori- 
zontal advection  of  u,v,T  and  h;  the  vertical  advection  of  u,v  and 
T;  the  Coriolis  acceleration;  the  centripetal  acceleration;  and 
the  pressure  gradient  force.  Next,  the  horizontal  and  vertical 
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diffusion  terms  are  calculated  using  u,v  and  T  at  time  level  n-1 
and  are  added  to  the  dynamical  terms.  The  coupling  procedure 
given  by  (6.1)  is  then  applied  to  these  accumulated  tendencies 
which  are  then  multiplied  by  2At  and  added  to  the  value  of  the 
variable  at  time  level  n-1.  This  partially  advances  the  fields  to 
the  new  time  level  n+1.  The  dynamic  stability  condition  (see 
above  and  in  Section  5)  is  then  imposed  and  the  variables  within 
the  mixed  layer  are  vertically  averaged  above  the  new  mixed  layer 
depth.  Finally,  the  tendencies  due  to  the  surface  fluxes  (A20) 
and  turbulent  mixing  (Section  4)  are  calculated  and  the  fields  are 
advanced  to  their  final  value  at  the  time  level  n+1. 

When  the  two  step  Euler-backward  integration  scheme  is  used, 
the  sequence  of  calculations  is  the  same  as  described  above  in 
both  the  "forward"  and  the  "backward"  steps.  The  only  difference 
is  that  the  diffusion  terms  are  no  longer  lagged  in  time  but  are 
calculated  from  the  same  fields  as  the  dynamical  terms. 
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TABLE  1 


Symbol       Value  Description 


*M 

2  x  108  cm2  s  * 

\ 

7   2-1 
4  x  10  cm  s 

KM 

.5  cm  s~ 

Km 

2  -1 
.5  cm  s 

Qt 

-4    -1 
2  x  10   deg 

To 

5°C 

po 

-3 
1.0276  gm  cm 

f 

6  x  10~5  s_1 

Ar 

4.5  km 

r 
m 

45  km 

R 

450  km 

H 

964  m 

At 

150  s 

Horizontal  eddy  viscosity 
Horizontal  eddy  conductivity 
Vertical  eddy  viscosity 
Vertical  eddy  conductivity 
Thermal  expansion  coefficient 
Reference  temperature 
Reference  density 
Coriolis  parameter 
Horizontal  grid  size 
Radius  of  maximum  winds 
Domain  extent 
Depth  of  ocean 
Time  step 
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FIGURE  CAPTIONS 

Fig.  1.   Forcing  functions  for  the  ocean  model,  each  normalized  by 
the  value  of  the  function  at  the  radius  of  maximum  winds. 
The  tangential  and  radial  stress  components  have  the  same 
form  (solid  line)  with  maximum  values  of  35.9  and  -12.9  dpa, 
respectively.  The  surface  heat  flux  (dashed)  has  a  maximum 
value  of  -840  Wn  . 

Fig.  2.   Schematic  illustration  of  the  embedded  mixed  layer  of  depth 
h,  located  within  the  k  th  level  in  the  circulation  model. 
See  Section  6  for  a  complete  discussion  of  the  embedding 
technique. 

Fig.  3.   Predicted  oceanic  response  after  12  h  to  the  hurricane 
forcing  in  Fig.  1  using  the  circulation  model  with  the 
embedded  mixed  layer.  Only  the  upper  150  m  is  shown,  with 
the  center  of  the  hurricane  (r=0)  on  the  left  and  each 
larger  tick  mark  on  the  abscissa  corresponding  to  45  km. 
(a)  Tangential  current  component  (cyclonic  flow,  dashed) 
with  contour  interval  of  0.1  m  s-1 .   (b)  Radial  current 
component  (inflow,  dashed)  with  contour  interval  0.1  m  s 
(c)  Vertical  velocity  (downward  flow,  dashed)  with  contour 
interval  of  1  m  h_1.  (d)  Temperature  contours  (°C)  . 

Fig.  4.   Similar  to  Fig.  3  except  at  24  h. 
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Fig.   5.       Temperature  deviations   (contour   interval  0.5°C,  negative 
values  dashed)    from  the  initial   state  predicted  by  the 
embedded  mixed  layer  model  after   (a)    12  h  and    (b)    24  h. 

Fig.  Al.     An  example  of  mixed  layer  deepening  and  cooling  by  entrap- 
ment.    The  original   structure  prior  to  entrainment   (solid 
line)    had  a  mixed  layer  depth  of  36  m,  mixed  layer  tempera- 
ture of  22.5  C  and  temperature  jump  at  the  base  of  the  layer 
of  approximately  2.6°C.     The  mixed  layer  deepened  by  entrain- 
ment to  below  40  m,   the  mixed  layer  temperature  decreased  to 
22.2  C  and  the  temperature  jump  at  the  base  of  the  layer 
decreased  to  approximately  2.4°C  (dashed  line)  .     Heat  is 
conserved,   so  area  A  =  area  B. 

Fig.  A2.     An  example  of  mixed  layer  shallowing.     The  original   tempera- 
ture structure  prior  to  shallowing    (solid  line)    is  the  same 
as  in  Fig.  Al.     In  this  example,   the  mixed  layer  shallowed 
to  a  new  depth  which  lies  above  the  k  th  level,  but  which 
will  not  be  apparent   in  the  temperature  structure   (dashed 
line)    until   there  is  a  downward  flux  of  heat  at  the  surface. 
The  new  temperature  structure  was  obtained  from   (A23)-(A25) 
and   is  such  that  area  C  +  area  D  =  area  E. 
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